clear all; clc;
% Hochrechnung auf eine Sekunde Bestrahlung, Soll:10^16  Photonen
simfaktor = 1e11; % 10^16 = simfaktor*(Anz. Photonen in Simulation)

% Von welcher Seite soll das Bild erstellt werden?
view = 'side'; % Auswahloptionen: top, front, side

%% Hinweis zu den Indizes:
% Die Indizes X, Y und Z stehen fuer die 3 Dimensionen in der 
%  Logdatei bzw. in der Simulation.
% Die Indizes 1 und 2 an den entsprechenden Variablen stehen fuer 
%  die 2 Dimensionen in den 2D Plots. Variablen mit X, Y und Z
%  sagen also nichts direkt ueber das Aussehen der Plots aus.

nDiv = 400;  % Anzahl Voxel in Logdatei und in Plot
nDivX = nDiv; nDivY = nDiv; nDivZ = nDiv;
% Laenge eines Voxels in mm
lenX = 400/nDivX; lenY = 400/nDivY; lenZ = 400/nDivZ;

tumorX = ceil(0.5*nDivX); 		% Position Mittelpunkt Tumor
tumorY = ceil(250/400*nDivY);
tumorZ = ceil(210/400*nDivZ);

labelX = 'Breite in mm';
labelY = 'Hoehe in mm';
labelZ = 'Tiefe in mm';

switch(view) % Abhaengig der Betrachtungsrichtung Werte setzen.
    case 'front'
        schnittebene = tumorZ;
        nDiv1 = nDivY; % Anzahl Voxel in beide Plotdimensionen
        nDiv2 = nDivX;
        len1 = lenY; len2 = lenX;
        tumor1 = tumorY; tumor2 = tumorX;
        label1 = labelY; label2 = labelX; 
    case 'top'
        schnittebene = tumorY;
        nDiv1 = nDivX; nDiv2 = nDivZ;
        len1 = lenX; len2 = lenZ;
        tumor1 = tumorX; tumor2 = tumorZ;
        label1 = labelX; label2 = labelZ; 
    case 'side'
        schnittebene = tumorX;
        nDiv1 = nDivY; nDiv2 = nDivZ;
        len1 = lenY; len2 = lenZ;
        tumor1 = tumorY; tumor2 = tumorZ;
        label1 = labelY; label2 = labelZ; 
end

% Tumor einzeichnen und weitere Werte fuer Visualisierung setzen
ppp = 100;	% Anzahl Pixel, die den Tumor markieren. 100 ist gut
tumorRadius = 20;	% mm
tumorSurf1 = zeros(2*ppp,1); 	% Arrays fuer Tumormarkierungen
tumorSurf2 = zeros(2*ppp,1);
for k = 1:ppp			% Tumormarkierungen mit Inhalt fuellen.
    tumorSurf1(k) = ceil(tumorRadius/len1*cos((k-1)/ppp*2*pi)+tumor1);
    helpT = (tumorRadius-len1)/len1*cos((k-1)/ppp*2*pi)+tumor1;
    tumorSurf1(k+ppp) = ceil(helpT);
    tumorSurf2(k) = ceil(tumorRadius/len2*sin((k-1)/ppp*2*pi)+tumor2);
    helpT = (tumorRadius-len2)/len2*sin((k-1)/ppp*2*pi)+tumor2;
    tumorSurf2(k+ppp) = ceil(helpT);
end

% Leere Sparcematrix zum Einlesen der Daten
matrix2D = sparse(nDiv1,nDiv2);

if(exist('filename', 'var') == 0) 	% Dateien einlesen
    [filename, pathname] = uigetfile('.log', 'MultiSelect', 'on');

    if ~iscell(filename)    % Wenn nur eine Datei gewaehlt wurde.
        filename = {filename};
    end
    
    kEnd = numel(filename); % Anzahl an einzulesender Dateien.
    
    for k=1:kEnd						% Dateien einzeln durchgehen.
        file = fopen(fullfile(pathname, filename{k}), 'r');
        input = fscanf(file,'%u\t%e\t%e\n',[3 inf]);
        fclose(file);
        
        % Anzahl Eintraege/"Zeilen" in aktuell eingelesener Datei
        entries = numel(input(1,:));
        
        for i=1:entries 		% Durch alle Eintraege durchgehen.
            % Jeder Eintrag entspricht einem Voxel und jedes Voxel 
            %  hat eine Nummer. Anhang dieser Nummer wird die 
            %  Koordiante des Voxels (nx, ny, nz) berechnet. Die 
            %  Voxelnummer ist in Spalte 1 des Arrays "input".
            % Es gilt: Voxelnummer = nz + ny*nDivZ + nx*nDivY*nDivZ
            nx = floor(input(1,i) / (nDivY*nDivZ));
            ny = floor((input(1,i) - nx*nDivY*nDivZ) / nDivZ);
            nz = input(1,i) - nx*nDivY*nDivZ - ny*nDivZ;
            
            switch(view) % Werte fuer Betrachperspektive setzen.
                case 'front'
                    n1 = ny+1; n2 = nx+1;
                    nSchnitt = nz+1; 		% Kopf wird invertiert.
                case 'top'
                    n1 = nx+1; n2 = nz+1; 
                    nSchnitt = ny+1;
                case 'side'
                    n1 = ny+1; n2 = nz+1;
                    nSchnitt = nx+1;
            end
            
            % Daten werden nur fuer die Schnittebene, die der 
            %  Betrachter spaeter als Plot sieht, rausgeschnitten.
            if nSchnitt == schnittebene
                matrix2D(n1, n2) = matrix2D(n1, n2) + input(3,i)*simfaktor;
            end
        end
        clear input; % clear memory
    end
end

fig = figure('Position', [0, 0, 800, 664]); % Figure erstellen
load cmap; colormap(cmap); % Colormap laden

set(gca, 'XLim', [0 nDiv2*len2]); set(gca, 'YLim', [0 nDiv1*len1]);

%clim = [1e-3 1e-1]; 	% anzuzeigendes Energieinterval Strahlung
clim = [1e1 1e3]; 		% anzuzeigendes Energieinterval Ionen

for k = 1:ppp 				% Tumor wird in Daten eingezeichnet.
    matrix2D(tumorSurf1(k), tumorSurf2(k)) = clim(2);
end

linear_axes = subplot(1,1,1); 	% Plotten
lin_plot = pcolor(linear_axes, matrix2D);

set(gca, 'CLim', clim);
colorbar('Yscale', 'log', 'FontSize', 14);

% Logarithmischer Plot (Matlab ausgetrickst)
set(linear_axes, 'Visible', 'off');
set(lin_plot, 'Visible', 'off');

log_axes = axes('Position', get(linear_axes, 'Position'));
log_plot = pcolor(log_axes, log10(matrix2D));
xlabel(label2, 'FontSize', 14);
ylabel(label1, 'FontSize', 14);
set(log_plot, 'EdgeColor','none');
set(gca, 'CLim', log10(clim));

axis square tight;

title('Radiation Dosimetry of Patient in Grey (Gy)', 'FontSize', 14);
